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Abstract. We present in this work results from numerical solutions, obtained by 
means of the direct simulation Monte Carlo (DSMC) method, of the Boltzmann and 
Boltzmann-Lorentz equations for an impurity immersed in a granular gas under planar 
Couette flow. The DSMC results are compared with the exact solution of a recent 
kinetic model for the same problem. The results confirm that, in steady states and 
over a wide range of parameter values, the state of the impurity is enslaved to that of 
the host gas: it follows the same flow velocity profile, its concentration (relative to that 
of the granular gas) is constant in the bulk region, and the impurity/gas temperature 
ratio is also constant. We determine also the rheological properties and nonlinear 
hydrodynamic transport coefficients for the impurity, finding a good semi-quantitative 
agreement between the DSMC results and the theoretical predictions. 

Keywords: granular matter, kinetic theory of gases and liquids, rheology and transport 
properties 
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1. Introduction 

The transport of granular matter has a growing interest for industrial and technological 
purposes. For this reason, it is convenient to study the behavior of simplified granular 
systems, both from an experimental and a theoretical point of view (see [l], for 
example, for a recent review on the field). Some of the phenomena of interest for 
industrial applications are the behavior of transport, diffusion, and segregation of 
grains depending on their different physical properties (mass, form, size, or inelasticity) 
[2l[3lia|5l|6l[71|8l[9l[T0l[ni[I3[T3l|ia[l5l[^ When the granular medium is dilute 
and vigorously shaken, the motion of grains resembles that of atoms or molecules 
in an ordinary gas and the near-instantaneous binary collisions prevail. Under these 
conditions, kinetic theory properly modified to account for the inelasticity of collisions 
provides a useful framework to analyze granular flows. It has been shown that the 
corresponding kinetic equation describing the statistics of this many-particle system 
can generate the so-called 'normal solution', in which all the spatial dependence occurs 
through the average fields [17]. In this situation, the fluxes can also be expressed 
as functions of the average fields and this results in a closed set of equations for the 
average fields that define a hydrodynamic description since it is formally equivalent to 
the hydrodynamic description in classical fluid mechanics [18] . Moreover, if the spatial 
gradients are supposed to be small enough, the hydrodynamics is Newtonian and the 
resulting equations are the Navier-Stokes (NS) ones [T71 [19]. Therefore, there have 
been attempts to determine the NS transport coefficients of granular mixtures in the 
low density regime [8l|20l[2T] and also at moderate densities [22| |23| [2^ |25| [26| l27t l28| 129] . 

On the other hand, the inelasticity in the collisions (i.e., the kinetic energy loss in 
interparticle collisions), introduces an inherent time scale that results, if a steady state 
needs to be maintained, in a minimum size of the gradients as a function of the degree 
of inelasticity [30]. This renders the NS approach not appropriate for most steady rapid 
granular flows. In fact, we already know that, for example, the simple shear flow for a 
granular gas [HI [32l |33] is inherently non- Newtonian [34] . In addition, previous works 
on the latter flow have shown that nonlinear effects can significantly modify segregation 
criteria with respect to NS hydrodynamic theories [15] . 

Obviously, the determination of non-Newtonian hydrodynamic profiles and 
transport coefficients is a prerequisite for a more complete description of transport and 
segregation of impurities immersed in a granular gas. This fact motivated a previous 
work by the authors, in which, via a kinetic model for a multicomponent granular gas 
[35] . we determined an exact analytical solution of the problem of the steady planar 
Couette fiow for an impurity in a low density granular gas [36]. This solution presents 
the advantage of including terms of all orders in the velocity and temperature spatial 
gradients, thus capturing nonlinear effects such as normal stress differences and a heat 
flux component normal to the thermal gradient. On the other hand, recent theory 
results on monocomponent granular gases have revealed new and interesting classes of 
flows that support the validity of a hydrodynamic description, without the restriction 
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Figure 1. Sketch of the planar Couette flow. The granular gas and the impurity 
are enclosed between two infinite parallel walls located at y = ±L/2, moving along 
the x-direction with velocities ±[//2, and kept at the temperature T^. In the steady 
state, the granular gas and the impurity move with the same flow velocity and show 
different temperature profiles, 12(2/) and respectively, but with constant ratio 

X = T1/T2. 

of small spatial gradients [STJ [381 ESI SOl El Hlj • It is thus interesting (a) to check if 
the hydro dynamic profiles found in the kinetic model description are shared by the true 
Boltzmann equation and (b) to gauge the degree of accuracy of the analytical solution 
derived in [36] . 

This has motivated the present work, where by means of the direct simulation 
Monte Carlo (DSMC) method [43] we obtain the numerical solution of the Boltzmann 
equation associated with a smooth hard sphere impurity immersed in a low density gas 
of inelastic smooth hard spheres under Couette flow. As in the case of the kinetic model, 
the DSMC solution is in principle not restricted to small spatial gradients |43] . 

Let us consider a set of identical smooth hard spheres (of diameter (T2 and mass 7712) 
that collide inelastically with each other. The inelasticity of collisions is characterized 
by a constant coefficient of normal restitution a2- We will consider that the number 
density of the system is sufficiently low that the typical contact times at collisions are 
much smaller than the typical time interval between collisions. In this low density 
regime, one can neglect the velocity correlations between the particles that are about to 
collide ('molecular chaos' hypothesis). In these conditions, as for an ordinary gas of hard 
spheres, we can describe statistically this granular gas through the Boltzmann equation 
conveniently adapted to take into account the inelasticity in the collisions [191 EH SH] • 
The granular gas is enclosed between two infinite parallel walls, here assumed to be 
both at the same temperature T^, moving with a relative velocity U. In this way, the 
sheared granular gas reaches a steady state with a non-zero velocity profile. The base 
laminar fiows in this geometry (see figure [1]) should be of the form U2 = U2^x{y) ^x- 
Moreover, a temperature profile T2{y) and a density profile n2{y) are present. We 
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introduce now another set of inelastic smooth hard spheres (of diameter ai and mass 
mi). The concentration of this new species is assumed to be neghgible, i.e., its density 
(ni) relative to that of the granular gas (77.2) tends to zero: xi = ni/n2 — (tracer 
limit). For this reason we call the 'impurity' granular species number 1. The inelasticity 
of a collision between a sphere of species 1 and a sphere of species 2 is characterized by 
a constant coefficient of normal restitution ai, which in general differs from a2- Due 
to the shearing from the boundaries, the impurity species will also reach a base steady 
flow of the same type as that of the granular gas, i.e., ui = ui^xiy) e^, Ti{y), and ni{y). 

From the exact solution to the kinetic model mentioned above, one finds that the 
hydrodynamic profiles satisfy the following properties. Regarding the excess component 
(host granular gas) |16] , (i) the pressure p2 = n2T2 is uniform, (ii) the (local) shear rate 
du2,x/dy scaled with respect to the (local) collision frequency 1^2 oc n2Tl^'^ is uniform, 
and (iii) the temperature T2 is a parabolic function of the fiow velocity M2,x- In the case 
of the tracer component (impurity), one finds that [36] (iv) the fiow velocity coincides 
with that of the gas, i.e., ui^x = ^2,a:, (v) the mole fraction ni/n2 is uniform, and (vi) 
the temperature ratio x = T1/T2 is also uniform. As we said, the main goal of the 
paper is to confirm these predictions by means of computer simulations. In addition, 
we will measure the momentum and heat fiuxes of the impurity to get the generalized 
transport coefficients and compare them with the analytical results derived from the 
kinetic model. 

The structure of the paper is the following. We formally describe the kinetic theory 
problem in section O where the generalized rheological and transport coefficients are also 
defined. In this section [2] we also briefiy recall the kinetic model results and describe the 
numerical method (DSMC). In section [3] we present the simulation data compared with 
the analytical solution of the kinetic model and discuss the results. Finally, we present 
the conclusions in section |H Additionally, we present in the appendix the theoretical 
hydrodynamic properties for the impurity from our previous work [36] . 

2. Theoretical description and numerical methods 

2.1. The Boltzmann description of the Couette flow for the granular gas and the 



We consider a granular gas composed of inelastic rf-dimensional hard spheres of diameter 
(T2, mass m2, and coefficient of normal restitution a2- In the low density regime, the 
corresponding velocity distribution function /2(r, v, t) (in the absence of gravity) obeys 
the Boltzmann equation 



where J22[v|/2,/2] is the (inelastic) colhsion operator. 

We add now the impurity particles of diameter cti and mass mi (species 1), which 
are present in a vanishing concentration. For this reason, we can assume that collisions 
between impurity particles themselves may be neglected and, in addition, the state of 



impurity 
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the granular gas (species 2) is not affected by the presence of impurities, so equation 
(12. ip is still valid. The velocity distribution function /i(r, v, t) for the impurity particles 
obeys the Boltzmann-Lorentz equation 

(5t + v-V)/i = Ji2[v|/i,/2], (2.2) 

where Jui^lfi, f2] is the corresponding (inelastic) collision operator, which is 
parameterized by the impurity-gas coefficient of normal restitution ai. Detailed 
expressions for Jjj[v|/j, fj] can be found in, for instance, [35] . 

The relevant hydrodynamic fields for both species are the number densities rij, flow 
velocities Uj, and granular temperatures Tj. They are defined by the relations 



n,{r,t) = J dv/,(r,v,t), (2.3) 
u,(r,t) = I dvv/,(r,v,t), (2.4) 

Ti{r,t) = ^^^^ j dv [v-u,(r,t)]V^(r,v,t). (2.5) 

In equation (12. 5p we have defined the partial temperatures Tj taking the velocities of 
species i relative to its mean value Uj. The usual choice, however, is to refer the velocities 
to the global mean flow velocity u |36]. Here, for convenience, we adopt the former 
choice. In any case, since u = U2 in the tracer limit, both choices are equivalent in the 
case of the granular gas. Additionally, the pressure tensor Pj and the heat flux for 
each species can be defined as 

Pj(r, t)=mi j dv [v - u,;(r, t)] [v - Ui(r, t)] /,(r, v, t), (2.6) 

q,(r,t) = ^ j dv [v-u,(r,t)]'[v-u,(r,t)]/;,(r,v,t). (2.7) 

Given that the microscopic state of the granular gas is independent of the 
microscopic state of the impurities, the mass, momentum, and energy balance equations 
for species 2 take the usual form [T7] 

Dtn^ + naV ■ U2 = 0, (2.8) 
Au2 + ■ P2 = 0, (2.9) 

AT2 + (V ■ q2 + P2 : Vu2) = -C2T2, (2.10) 
dn2 

where = (9j + U2 ■ V is the material time derivative and 

^21, 2 



C2 = -^ / dvi;V22[v|/2,/2] (2.11) 



is the cooling rate of the granular gas. 
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Now we assume that the system (granular gas plus impurities) is subject to the 
planar Couette flow. For steady states {dt = 0), and given the geometry of the problem 
{dx = dz = 0), we obtain 

dyP2,xy = dyP2^yy = 0, (2-12) 

dyq2,y + P2,xydyU2,x = -^C2'^2T2. (2.13) 

Thus far, all the equations in this section are formally exact in the framework of 
the Boltzmann equation. Based on previous results [36] |46] derived from the kinetic 
model described below, we expect that the hydrodynamic flelds of the gas in the bulk 
domain of the system have the forms 

P2 = ^2X2 = const, (2.14) 

— dyU2^x = a = const, (2-15) 
1^2 

1 / 1 ^ ^ 



2m2Pr \ z/2 



dy] T2 = -7 = const. (2.16) 



Here, Pr = {d — l)/d is the conventional Prandtl number [18] and z/2 n2Tl^'^ is a 
characteristic collision frequency. For the sake of concreteness, henceforth we will take 

o_-(d-l)/2 

In equation fl2.15p . the constant a represents a dimensionless shear rate. It plays 
the role of a Knudsen number associated with the velocity gradient. The constant 7 
in equation fl2.16p is a dimensionless parameter (henceforth called thermal curvature 
coefficient) characterizing the curvature of the temperature profile as a consequence of 
both the viscous heating and the coUisional cooling. As a consequence, 7 must depend 
on both the shear rate a and the coefficient of restitution a2- It is interesting to note 
that, from equations f l2.15p and f l2.16p . one finds 

r2(y) = r2(0)-Pr^<(y), (2.18) 



a 



where we have considered that M2,x(0) = and the temperature profile is symmetric. 
Equation fl2.18p implies that, if eliminating y between T2 and M2,x5 the temperature is a 
linear function of u\^. 

With respect to the state of the impurities, we assume that it is enslaved to that 
of the gas. This means that 



Ul = U2, (2.19) 

Xi = — = const, (2.20) 

n2 

X = ^ = const. (2.21) 
J 2 
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P,., = -v:(a)^^. (2.22) 
e.Ja) = 1^, 9,Ja) = (2.23) 



The hypothesis (12.191) implies that there is no diffusion of the impurities with respect to 
the gas particles. Equations (I2.20p and (I2.2ip . together with equation (I2.14p . imply that 
Pi = riiTi = const. In summary, the hydrodynamic profiles of the system are provided 
by equations ^J^-^J6^ and (ICTjl -f iOTj) . 

In order to characterize the non- Newtonian properties, it is convenient to introduce 
the following five dimensionless rheological factors [361 SIl 1121 06] 

)y 

P P 

d + 2 UiTidTi 

d + 2 UiTidTi 

where 1/2 is defined by equation (I2.17P and 

^^- id + 2nd/2r ^'^\^^^) • (2.26) 
Note that, in the case i = 1, the definitions of rj*, A*, and (p* in equations (I2.22p . 
(I2.24p . and (I2.25P slightly differ from those in [36]. The nondimensionalization of the 
generalized shear viscosity rj* and thermal conductivity A* is such that 773 = Ag = 1 
for an elastic gas (02 = 1) in the NS regime (a — 0) [18]. Note that, while 77* and 
A* are generalizations of NS transport coefficients, the functions Oi^^, 9i^y, and 0* are 
generalizations of Burnett transport coefficients [HI HTl HE] . 

Of course, if the impurities are mechanically equivalent to the gas particles (i.e., 
mi = 7712, CTi = (T2, and ai = 02), the transport coefficients associated with both species 
coincide. 

Taking into account the constitutive forms (I2.22p and (I2.24p for the granular gas 
{i = 2), as well as equations (I2.14p - (l2.16p . the exact balance equation (I2.13P becomes 

^;a'-id + 2)x;^ = ^C2, Q^-- (2.27) 

Note that in the elastic case {Q = 0) and in the NS limit (i.e., 773 — )■ 1 and Ag — ?■ 1) one 
has 7 = a'^/{d + 2) [39]. In general, 7 depends on both a and 0^2, its sign depending 
on the competition between viscous heating and inelastic cooling. If viscous heating 
dominates (i.e., r/gO^ > dQ/2), then 7 > 0. On the other hand, 7 < in the opposite 
situation (i.e., r/ja^ < dQf^)- Both effects cancel each other (and thus 7 = 0) at a 
threshold shear rate ath given by 

= t^Sr-y (2.28) 

Since 7 = at a = ath, it follows from equations fl2.16p and fl2.2ip that u~^dyTi = const. 
This condition, together with TZjTj = const and equations (12.240 and (12.250 . implies that 
the heat flux is uniform at a = ath S21 ED] both for the impurity and the host gas. 
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2.2. The kinetic model description of the Couette flow for the granular gas and the 
impurity 

The mathematical complexity of the Boltzmann and Boltzmann-Lorentz equations (12.11) 
and (12. 2p prevents one from obtaining exact solutions. This has motivated the proposal 
of simpler kinetic models, most of them inspired on the well-known Bhatnagar-Gross- 
Krook (BGK) kinetic model for ordinary gases [51]. Here we consider the following 
BGK-type kinetic model [351 [52] : 



1 + a- 

^»2[v|/i, /2] -kd - V , (/i 



[(v - Ui) fi^ 



where 



/.2(V) 



rii 



rrii 



2nT, 



i2 



d/2 



exp 



2T, 



Ui2 



i2 



(2.29) 



(2.30) 



is a reference distribution function. In the case of the granular gas (z = 2), T22 = T2 
and U22 = U2, so /22 is the local equilibrium distribution function. In the case of the 
impurity particles {i = I) [35] . 



12 



2/i 



U12 



T2-ri 



yUUl + U2 



[Ui - U2j 

2d 



m2 + 



Ti 



Ti/mi + T2/m2 



where 



mi 
/i = — 
m2 



is the mass ratio. In equation (12.291) 
d + 2 



C2 



Ad 



'1 



«2)'^2 



is the cooling rate (12. lip evaluated in the local equilibrium approximation, while 



Ci 



d + 2 



2d 



/i)2 



miT2 ^ 3 mi 



Ui - U2 



(2.31) 
(2.32) 

(2.33) 

(2.34) 

(2.35) 



m2Ti ' 2d Ti 

is the impurity cooling rate. Finally, the factor kd can be chosen to optimize agreement 
with the Boltzmann description. In particular, the choices kd = 1 and kd = Pt = 
[d — l)/d reproduce the NS shear viscosity and thermal conductivity coefficients, 
respectively, of the gas in the elastic limit [TH [51] . A third criterion to fix the factor kd 
is to require that the coUisional momentum transfer of the impurities be the same for 
the kinetic model as for the true Boltzmann-Lorentz equation [351 EHj- This results in 
kd={d + 2)/d. 

In a recent work [36], we solved the kinetic model equations for the granular gas 
and impurity in the steady Couette fiow. The resulting profiles agree with the forms 
(l27[4ll -( [2:T6ll and ([CTD - dOTl) . Moreover, the solution gives 7, x, Vi, di,y, AJ, and 
(pl as functions of a, ai, a2, fi, and u = (Ti/(T2- Their explicit forms are displayed in the 
appendix. 
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Figure 2. This figure illustrates the verification in the bulk region of the three 
hypotheses on the stationary Couette flow for the granular gas: (a) the hydrostatic 
pressure p2 is constant, (b) the local shear rate a is also constant, and (c) the 
temperature T2 is a linear function oi ^. Two values of the coefficient of restitution 
are considered: 02 = 0.9 (Q, with L = 23.23 and a = 0.443) and a2 = 0.8 (□, with 
L = 15.48 and a = 0.571). Lines in (c) stand for linear fits to DSMC data. 



2.3. Numerical methods (Monte Carlo simulations) 

As said in section [H we have solved numerically the Boltzmann and Boltzmann-Lorentz 
kinetic equations for the hard-sphere {d = 3) granular gas and impurity [equations (12 .ip 
and (12.2p . respectively] by means of the DSMC method. Originally devised for elastic 
gases in the low density limit [13], this method has been successfully implemented also 
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Figure 3. This figure illustrates the verification in the bulk region of the three 
hypotheses on the stationary Couette flow for the impurity: (a) no mutual diffusion 
exists, i.e., ui^xiy) — U2,x{y), (b) the mole fraction ni/n2 is constant, and (c) the 
temperature ratio x — T1/T2 is also constant. Two cases are considered: ai ~ a2 = 
0.9, /.t = 2, w = 1 (O, with L = 23.23 and a = 0.443) and m = a2 = 0.8, = 0.5, 
Lo — 1 (□, with L — 15.48 and a ~ 0.571). In panel (a) the crosses and triangles 
correspond to the granular gas (i = 2), and the symbols corresponding to a = 0.8 have 
been lifted to avoid overlap with the symbols corresponding to a = 0.9. In panel (b) 
the density rii is normalized with respect to its spatial average value Hi. 
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for inelastic hard spheres [H] and moderately dense gases (Enskog kinetic equation 
[531 El [55]). 

The implementation of the algorithm for granular gases has been described in more 
detail elsewhere (see, for instance, [12]). We will recall only that it consists of small 
time steps (in the scale of the characteristic mean free time), each one having two 
basic stages: (i) free streaming and (ii) stochastic interparticle (binary) collisions. The 
system is divided into small cells (in the scale of the characteristic mean free path). 
In the DSMC method, the number Ni of simulated particles is a statistical technical 
parameter that, in contrast to the molecular dynamics case, does not need to coincide 
with the actual number of physical particles in the system. We have used the same 
number of simulated particles (A^i = A''2 = 2 x 10^) for both species. On the other hand, 
since the physical impurity concentration is assumed to be negligibly small, only 2-2 
collisions are considered in the evolution of the granular gas and only 1-2 collisions are 
considered in the evolution of the impurity particles. Therefore, after a collision of type 
2-2, the velocities of both particles are changed, while only the velocity of the impurity 
particle is changed after a collision of type 1-2. 

As in a previous work [42], we perform two averages for steady states: (i) a first 
spatial average over neighbor simulation cells, taking care that this coarse-grained cell 
size is not larger than the typical length over which hydrodynamic fields vary, and (ii) 
a time average for each coarse-grained cell since the microstates of the simulation are 
stored iteratively many times during the stationary state in the simulation. 

Since the parameter space (mass ratio /i = mi/m2, size ratio oj = oxjo^-, coefficients 
of restitution a\ and a^-, and shear rate a) is quite large, we have focused on a few 
representative cases. As in our previous work on the BGK-type model [36], we have 
analyzed cases with a common coefficient of restitution (ai = 0^2 = a) and equal 
diameter [bj = 1) for d = 3 (spheres). Three different values of the mass ratio 
(/i = 2, 1, 0.5) and of the coefficient of restitution [a = 1, 0.9, 0.8) have been considered. 
Thus, the elastic limit (a = 1) as well as the properties of the granular gas (/i = 1) are 
included as particular cases. For each one of the nine combinations of the pair (/i, a), 
we have made series of simulations in shear rate a by keeping fixed the wall velocity 
difference U = 10, but varying the wall distance in the range L = 2.5-30. We will 
use hereafter nondimensionalized quantities with the following choice of units: m2 = 1, 
T2(±L/2) = 1, and U2{±L/2) = 1. Also, we define the density levels Ui = 1, where the 
bar denotes a spatial average across the system. 

In section |3] we compare the DSMC numerical solution of the Boltzmann description 
with the analytical solution of the BGK-type kinetic model. But, before that, it is 
appropriate to show DSMC data confirming that the hypotheses on which our theoretical 
description relies on [equations f l2.14p . f l2.15p . f l2.18p -( l2.2ip ] are indeed valid. As an 
illustration, figures [2] and [3] display the DSMC profiles for the cases {fj,, a) = (2, 0.9) 
and (/i, a) = (0.5,0.8). From figure |2] we observe that the hydrostatic pressure p2 of 
the granular gas, except for small infiections near the boundary layers, is fiat, the local 
shear rate a is indeed constant throughout the system, and the temperature T2 is a linear 




a 



Figure 4. Threshold values ath for the shear rate a vs the coefficient of normal 
restitution for spheres (d = 3). The lines represent the theoretical results given by 
equation (|A.10p with three difFerent choices for the parameter kd- kd = 1 (solid line), 
kd = Pt = {d—l)/d (dashed line), and kd = {d + 2)/d (dotted line). Symbols stand for 
DSMC simulation data with AT = (■, this work), AT = 2 (x, glllig), AT = 10 
(+, [311112]), and simple shear flow (Q. [56]). 

function of u^^. Figure [3] confirms the remaining hypotheses fl2.19p - fl2.2ip . specific for 
the granular impurity. The "enslaving" condition Ui = U2 is fulfilled with a very high 
degree of accuracy for all points in the system, even in the boundary layers. Moreover, 
the relative impurity concentration and the temperature ratio are practically constant. 

3. Results and discussion 

3.1. The threshold shear rate 

We first consider the threshold value of the shear rate ath at which the thermal curvature 
parameter 7 vanishes. As discussed below equation f l2.27p . the value a = ath is especially 
important since it corresponds to an exact balance between viscous heating and inelastic 
cooling, giving rise to = const. In the geometry sketched in figure [H where both 
walls are maintained at the same temperature, the value 7 = implies a constant 
temperature T2 and thus the Couette flow becomes equivalent to the well-known simple 
shear flow |3T1 |32l |33]. More in general, when the walls are allowed to have different 
temperatures (AT = T^+fT^j^ — 1 7^ 0), the limit case 7 = defines, in the parameter 
space {a, a, AT}, a surface that has been shown recently [411 142] to represent a special 
class of granular flows, including the conventional Fourier flow of an elastic gas. This 
generalized class of flows (that we called "LTu" because the temperature is a linear 
function of the flow velocity) can be theoretically described in a single hydrodynamic 
theory frame, both for elastic and granular gases. 

Figure H] shows the a-dependence of the threshold shear rate ath, as predicted by 
the BGK-type kinetic model, equation flA.lOp . with the three choices of mentioned 
before, namely = 1, k^ = {d — l)/d, and kd = {d + 2)/d. The simulation data obtained 
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Figure 5. Thermal curvature coefficient 7 as a function of the square of the shear rate, 
a^. From now on, in all figures, lines stand for the BGK-type theory (with fc^^ = 1) 
and symbols stand for DSMC data. Three series are plotted: a — 1 (solid line and 
circles), a = 0.9 (dashed line and squares), and a = 0.8 (dotted line and triangles). 
The threshold levels a^j^ for the shear rate in the cases a = 0.9 and a = 0.8 are marked 
with dotted-dashed vertical lines. 



here for a = 0.8 and 0.9 with AT = 0, as well as those of |12] for a wider range with 
AT 7^ 0, are also included in figure IH We clearly observe that the best agreement is 
achieved with the choice = 1. Although we have checked that either kd = {d — l)/d 
or kd = {d + 2)/d may provide a better agreement with simulation for some of the 
other quantities, henceforth we will adopt the choice /c^ = 1 as a convenient compromise 
between simplicity and accuracy. 

In the rest of this section, we will compare DSMC results from the Boltzmann 
equation with the analytical results derived from our BGK-type model, representing 
the relevant hydrodynamic properties as functions of the shear rate (for a > Oth), so we 
can analyze to what extent the nonlinear theoretical description is reliable, at least at 
a qualitative level, for these very far away from equilibrium steady states. 

3.2. Thermal curvature coefficient 7 and temperature ratio x 

In figure [5] we can see that the theoretical thermal curvature parameter shows a good 
agreement with DSMC data. The agreement is quantitatively very good near the 
threshold Oth, but the theory tends to overestimate 7 as the shear rate increases. In any 
case, we can observe in figure [5] the reliability of our non-Newtonian hydrodynamic 
description in the context of the BGK-type kinetic model. As usual, the thermal 
curvature parameter decreases for decreasing shear rate, until it reaches a threshold 
value 7 = at which we obtain the simple shear flow [33] or, more generally, the LTu 
class flow H2| 150] . We should recall that the BGK solution is not mathematically 
well defined for states with 7 < (see, however, [57] for an analytical continuation). In 
principle, it should be also possible to find nonlinear Couette flows from the Boltzmann 
equation (it has already been shown that these states do exist in the quasielastic limit 
[30]). This region would correspond to points below the LTu surface [HI H2] . 
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Figure 6. Temperature ratio x = T1/T2 as a function of the square of the shear rate, 
a^. Two mass ratios are considered for each value of a: /i = mi/m2 = 2 (dashed hues 
and circles) and /.t = 1/2 (dotted lines and squares). The threshold levels a^j^ for the 
shear rate in the cases a = 0.9 and a = 0.8 are marked with dotted-dashed vertical 
lines. 

In figure [6] we plot the results for the temperature ratio. We may confirm that the 
Boltzmann equation solution follows the same general trends as the BGK-type model 
[36] . In particular, the impurity has a higher (lower) temperature than the granular 
gas if its mass is larger (smaller) than that of a gas particle, this effect being more 
pronounced as the shear rate increases. On the other hand, the BGK-type model tends 
to exaggerate these effects. 

3.3. Generalized transport coefficients 

In figure[7]we present the results for the shear viscosity 77^ As we can observe, the kinetic 
model predicts that tjI is practically independent of the mass ratio /i. This feature is also 
shared by the DSMC data, except near Oth, where the values of rjl increase with the mass 
ratio. Apart from this, the kinetic model not only successfully captures the decrease 
of the shear viscosity with increasing shear rate (shear thinning), but also exhibits a 
generally good quantitative agreement. 

The normal stress coefficients 9i^x and 9i^y are displayed in figures [8] and M, 
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Figure 7. Reduced shear viscosity 77* as a function of the square of the shear rate, a^. 
Three mass ratios are considered for each value of a: fi = 7711/7712 = 2 (dashed hnes 
and circles), fJ, = 1 (solid lines and triangles), and n = 1/2 (dotted lines and squares). 
The threshold levels a^j^ for the shear rate in the cases a = 0.9 and a = 0.8 are marked 
with dotted-dashed vertical lines. 

respectively. While 6i^r^ > 1, we observe that 6i y < 1, this anisotropic effect increasing 
with increasing shear rate, with increasing mass ratio, and with increasing coUisional 
dissipation. The kinetic model correctly accounts for these trends. At a quantitative 
level, the agreement is fairly good for /i < 1 (especially in the case of Oi^x)- However, 
in the case of heavy impurities (yU > 1), the kinetic model clearly underestimates the 
deviations of 6i^x and 6i^y from unity. 

The heat flux transport coefficients XI and (pl are shown in figures [ID] and [HI 
respectively. The theoretical curves for A^^ and, especially, for (pl are more sensitive to 
the value of the mass ratio n than those for rjl, showing that the heat fiux transport 
coefficients increase with increasing fi. These features are generally confirmed by the 
DSMC data for (pl but in the case of the infiuence of fi is less clear, except near 
the threshold shear rate. Apart from that, the kinetic model correctly describes the 
decrease of AJ with increasing shear rate as well as the change in the dependence of (pl 
on a as one goes from the elastic case (a = 1) to the inelastic ones {a = 0.9 and 0.8). It 
is interesting to remark that the cross thermal conductivity coefficient (pl (an obvious 
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Figure 8. Normal stress cocfEcient di ^ as a function of the square of the shear rate, 
a^. Three mass ratios are considered for each value of a: /i = nii/ni^ — 2 (dashed lines 
and circles), = 1 (solid lines and triangles), and ^ = 1/2 (dotted lines and squares). 
The threshold levels a^jj for the shear rate in the cases a = 0.9 and a = 0.8 are marked 
with dotted-dashed vertical lines. 

non-Newtonian effect) can become larger than the generahzed NS thermal conductivity 
coefficient XI. This effect was already observed in the special case of LTu ffows in one- 
component systems H2] . Now, both the kinetic model and the DSMC data show 
that the inequality (pl > XI (i.e., \qi,x\ > \(li,y\) becomes more pronounced as the shear 
rate or the mass ratio increase. 

4. Conclusions 

We have analyzed in this paper the properties of the Couette ffow for a granular impurity 
immersed in a granular gas. We have focused on the region of high shear rates, i.e., 
a > Oth, where viscous heating prevails over inelastic cooling and thus the thermal 
curvature coefficient 7 defined by equation f l2.16p is positive. This corresponds to states 
in the region above the LTu surface described in a previous work [HI |12], where an 
analytical solution from a BGK-type model is available [36] and, more importantly, 
nonlinear effects are dominant. We have confirmed that the hypotheses made in order 
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Figure 9. Normal stress coefficient 9i^y as a function of the square of the shear rate, 
a^. Three mass ratios are considered for each value of a: ^ = mi/m2 — 2 (dashed lines 
and circles), A* = 1 (solid lines and triangles), and = 1/2 (dotted lines and squares). 
The threshold levels a^-^ for the shear rate in the cases a = 0.9 and a = 0.8 are marked 
with dotted-dashed vertical lines. 

to obtain the BGK-type theoretical solution in a previous work |36l |16] also apply for 
the numerical solution of the inelastic Boltzmann and Boltzmann-Lorentz equations: (i) 
the pressure p2 of the granular gas is uniform, (ii) the local shear rate du2^xldy divided 

1 /2 

by the local collision frequency z/2 oc n2T2 is uniform, (iii) the temperature T2 of the 
granular gas is a quadratic function of the flow velocity ^2,^, (iv) the flow velocity of the 
impurity coincides with that of the gas, i.e., Ui^^ = ^2,2:; (v) the mole fraction ni/n2 of 
the impurity is uniform, and (vi) the impurity/gas temperature ratio x = T1/T2 is also 
uniform. The fourth hypothesis on the absence of mutual diffusion is fulfilled with an 
especially high degree of accuracy, since we were unable to measure any non-negligible 
difference between the profiles ui^^^y) and U2^x{y) over a wide range of situations [see 
figure m^a) for reference] . This is important since it delimits well the properties of the 
Couette flow for the impurity, which has important consequences in applications such 
as segregation [15]. 

Apart from the validation of the hydrodynamic hypotheses (i)-(vi), the DSMC 
results have shown a good semi- quantitative validity of the theoretical predictions as 
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Figure 10. Reduced thermal conductivity A* as a function of the square of the shear 
rate, a?. Three mass ratios are considered for each value of a: /i = mi/ 1712 = 2 
(dashed lines and circles), /i = 1 (solid lines and triangles), and ^ = 1/2 (dotted lines 
and squares). The threshold levels a^i^ for the shear rate in the cases a = 0.9 and 
a = 0.8 are marked with dotted-dashed vertical lines. 



regards the dependence of non-Newtonian properties (temperature ratio, normal stress 
differences, generalized shear viscosity, and generalized thermal conductivities) on the 
shear rate, the mass ratio, and the coefficient of restitution. 

Although we have limited our study to the region a > ath (i-e., 7 > 0), preliminary 
results seem to indicate that the special LTu state (a = ath or 7 = 0) also exists for 
the granular impurity (and exactly at the same point that occurs for the granular gas). 
This, together with the fact that the impurity in any case fulfills also the hypotheses 
(i)-(iii) for the granular gas, suggests that we can perform a classification of Couette 
flows for the impurity analogous to that for the granular gas [30l HI]. Therefore, an 
interesting way in which our work can be extended consists of studying Couette flows 
for non-equal temperature walls and in the region below Oth- Theoretical and numerical 
work on this subject is ongoing. 
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Figure 11. Reduced cross thermal conductivity (pl as a function of the square of the 
shear rate, a^. Three mass ratios are considered for each value of a: /i = mi/m2 = 2 
(dashed lines and circles), fJ. = 1 (solid lines and triangles), and /i = 1/2 (dotted lines 
and squares). The threshold levels a^j^ for the shear rate in the cases a = 0.9 and 
a = 0.8 are marked with dotted-dashed vertical lines. 
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Appendix A. Transport properties of the impurity from the BGK-type 
kinetic model 

We present in this appendix the explicit expressions for the relevant hydrodynamic 
properties, as extracted from a previous work [36] . Those expressions are given in terms 
of some mathematical functions that we define below. 
First, we introduce the functions 

POO 

Fo,m{y, = / dw e-(i+^)"'w;™Xo {Q{w, y, z)) , (A.l) 

^0 
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1 

2 Jo 



e{w,y,z) 



(A.2) 



where 



Xo(e) = v^ee®'erfc (9) - 1, Xi(0) = 6^ + 2e2)e®'erfc (6) - 20j , (A.3) 

1 z 



Q{w,y,z) 



(A.4) 



2V2yi-e 

In equation (]A.3p . erfc(x) is the complementary error function. Note that, because of 
the square root in equation (lA.4p . the functions FQ^m{y,z) and Fi^rn{y,z) are only well 
defined for y > 0. Next, we define 

'd+1 



poo 

Jo 



■-X,iQ{w,y,z)) + YiQiw,y,z)) 



H{y.z) 



dwe-^^+'i'^'"w'^Y{Q{w,y,z)), 



where 



r(e) = 2(1 + e^) - v^e(3 + 2e2)e« erfc (6) 



, (A.5) 
(A.6) 

(A.7) 



Now that we have introduced the above functions, let us display the expressions for 
the solution of the kinetic model. First, the thermal curvature coefficient 7 is given as 
a function of the reduced shear rate a and the coefficient of restitution 0:2 through the 
implicit equation: 



d 



C2C2 2c2a^ 



1 + C2Q (1 + C2Qf 



2^1,0(0^7, C2C2) + ^^^"0,0(0^, C2C2) 
+ cla^ [2Fi,2(c^7, C2C2*) + ^0,2 (ch, C2C2 



where 



Ci 



(A.8) 



(A.9) 



Since the functions ( lA.ll) and (I A. 2 1) are not defined for negative the representation 
(lA.SP exists only for 7 > or, equivalently, for a > Oth, where the threshold value ath 
of the shear rate (corresponding to 7 = 0) is [361 HH] 

'^-C*(1 + C2C2)'- (A.IO) 



2 



2Co 



In the case a = ath the viscous heating is exactly balanced by coUisional cooling and 
the heat flux becomes uniform H2| 150] . 

Next, the temperature ratio x = T1/T2 is obtained from the implicit equation 



d 



T12 1 + Ci 



2a' 



:i+ci)= 



2Fi,o(7,Ci) + rfi^o,o(7,Ci^ 



2Fi,2(7,Ci) + i^o,2(7,Ci) • (A.ll) 
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Here, 



a = — Cifl, 
1^1 



^ _ulTi2X 2 
7 = — Ci7, 
v( Ti /i 



~ ciCi d + 2 /i + x ,^ 2^ 



(A.12) 



(A.13) 



1^1 2d + 

In equation (1A.13I) use has been made of equation (]2.35p with ui = U2. Moreover, from 
equations (l^TTD . flT^ . and f lCT]) . we have 



1^1 



2/i 



Ti2 2/i(l-x) 



.1 + 

We recall that u; = oxlo2- 

Finally, the transport coefficients defined by equations fl2.22p - p.25p are 



(A.14) 



^1 = ci 



^12 



1 + Ci)^ 

d- 1 
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t7l,j; - — 



Ti Li + Ci 
1 T121 



Fo,i(7,Ci) + 2Fi,i(7,Ci 
(rf-l)Fo,o(7,Ci) + 2Fi,o(7,Ci) 

1 + Ci 

Fo,o(7,Ci) + 2Fi,o(7,Ci 



+ 2 Ti 7 



*~2 

?7ia — ci; 



J 1 



2ci 



ci + 2 V Ti 



12 



^27 



G(7,Ci) + a'i^(7,Ci) 



(A.15) 
(A.16) 
(A.17) 
(A.18) 
(A.19) 



where we have taken into account that Pr = 1 in the BGK model [51] . 

The transport coefficients for the granular gas are easily deduced from equations 
( IA.15p - (lA.19p by the formal replacements 1 — >■ 2, x — 1, /i — 1, and w — )■ 1. It is then 
easy to check that equation flA.lSP reduces to equation fl2.27p . 

To conclude this appendix, let us write the results in the limit a — ath, i-e., 7 — 0. 
The temperature ratio is given by the physical root of the quartic equation [36] 

1 \ 2Hi , 

(A.20) 



d{^ 



12 1 + Ci 



0, 



:i+Ci 



where ath is obtained from equations (lA.lOp and (IA.12p . Once x is known, the transport 
coefficients are [36] 

T12 1 



m = ci- 



Ti (1+Ci)^ 

9l,, = d-{d-l)9l,y, 



T12 1 



^1 1 + Ci 



(A.21) 
(A.22) 
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Tj 2 + 7Ci + 6Cf 



1 + 



6 12 + 42Ci + 37Ci2^2 
— flti 

(2 + 7Ci + 6Cf)2 



(A.23) 



^ (Til 
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2 



4 + 7Ci 



8 + 28Ci + 25Cf 
(2 + 7Ci + 6g)2 



+ 2 V Ti 



(2 + 7Ci + 6Cf )2 



ath c? + 4 + 18 




(A.24) 



Note that equations (lA.lOp and (lA.2ip (when the latter is particularized to the case 
of an impurity mechanically equivalent to the particles of the gas) are consistent with 
equation f lZ^ . 
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